library(mgcv); library(data.table); library(gratia); library(marginaleffects); library(ggplot2); library(dplyr); library(tidyr); library(zoo); library(slider); library(here);library(broom); library(readr);library(lme4)
data=read_tsv("G:/My Drive/cesab bioforest/Data/data.tsv")
sites=c("Corinto" ,"Jari" , "Lesong" ,"Mbaiki","Misiones", "Paracou",
"SUAS" ,"Sungai Lalang", "Ulu Muda")
setDT(data)
unique_plots <- data[site %in% sites, .(plots = unique(plot)), by = site]
setnames(data, "plot", "subplot")
data[, plot := gsub("_.*$", "", subplot)]
load("../clim_by_site_total.rda")
plots_df <- do.call(rbind, lapply(names(clim_by_site), function(site_name) {
data.frame(
site = site_name,
plot = unlist(lapply(clim_by_site[[site_name]], function(block) block$plot)),
stringsAsFactors = FALSE
)
}))
census_anom_all_sites<- list()
all_sites_results <- list()
for (s in sites) {
cat("Processing", s, "\n")
## ---- subset obs
obs_site <- data[site == s]
#obs_site2=obs_site[!duplicated(obs_site)]
obs_site2 <- obs_site[rel_year >= 0]
obs_site_wide <- dcast(
obs_site2,
treatment+site + subplot + plot+ year ~ variable,
value.var = "value"
)
## ---- climate stack
clim_all <- rbindlist(
lapply(clim_by_site[[s]], function(x) x$climate[[1]]),
idcol = "plot"
)
## ---- baseline
baseline <- clim_all[year >= 1981 & year <= 2010]
clim_stats <- baseline[, .(
mean_tmax = mean(tmax),
sd_tmax = sd(tmax),
mean_vpd = mean(vpd),
sd_vpd = sd(vpd),
mean_srad = mean(srad),
sd_srad = sd(srad)
), by = month]
clim_all <- merge(clim_all, clim_stats, by = "month", all.x = TRUE)
## ---- standardize
clim_all[, `:=`(
z_tmax = (tmax - mean_tmax) / sd_tmax,
z_vpd = (vpd - mean_vpd) / sd_vpd,
z_srad = (srad - mean_srad) / sd_srad
)]
## ---- match census years - needs to be by plot, as many sites dont have census for every plot every year
inv_years <- obs_site_wide[, .(year = sort(unique(year))), by = plot]
clim_all[
inv_years,
year_bin := i.year,
on = .(plot, year),
roll = TRUE
]
## ---- aggregate per census
census_anom <- clim_all[, .(
mean_z_tmax = mean(z_tmax, na.rm = TRUE),
mean_z_vpd = mean(z_vpd, na.rm = TRUE),
mean_z_srad = mean(z_srad, na.rm = TRUE)
), by = .(plot, year_bin)]#Because I group by year_bin, all pre-first-census months collapse into NA
## ---- merge obs + climate
merged_site <- merge(
obs_site_wide,
census_anom,
by.x = c("year","plot"),
by.y = c("year_bin","plot"),
all.x = TRUE
)
## ---- store
census_anom_all_sites[[s]] <- census_anom
all_sites_results[[s]] <- merged_site
}